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ABSTRACT 


The behavior of energetic particles in the solar system is 
described by a well known Fokker-Planck equation. Although 
analytic methods yield insight into the nature of its 

Jw V e T C * a A ly in the diffusion regime, calculations 
that go beyond diffusion are very complicated. Under these 
circumstances, numerical computations offer the only feasible 
way to obtain concrete results. However, the reliability of 
these calculations is of concern, because numerical methods 
are notorious for their errors and artifacts. To address 

i!n 8 n C n 0 r! Ce rt’ ^ Mlne P roblem of classical 
transport theory has been analyzed with the aid of three 

different numerical methods. These are: 


1. ^e method of eigenfunctions in which the distribution 
“J l0 J is approximated by a sum of eigenfunctions of the 

3 lng operator. Its complexity limits the practicality 
this approach, but it is closely related to the analytic 
and classical approaches. J 

2. Numerical solutions of a finite-difference equation. 
^b?le errors? 0811 praCtical approach, but it is subject to 

D * rect simulation of the scattering and streaming of 
individual particles with the aid of Monte Carlo methods. 

The accuracy of this approach is limited by statistical 
onsiderations , but it is closely related to the physics. 

If proper precautions are taken to ensure its validity, the 
second method gives results that are in precise quantitative 
agreement with those of the first. Results of Monte ?ar!o 

effe^! ti0 w a f 6 T accurate enough to define non-dif fusive 

»2 e d“f ft SL*SJ^. d “ , “ 1 “ re8lm<1 - they 

and Wigner (1958), the Milne 
J , ! the , “f ha cone of classical transport theory". It Is in 

» o r ° e a “1 a Ce8t f ° r ne ” '‘eialap-ants that this famous 

dfffislon of particles M ° re s Pacifically, the steady-state 

a fusion of particles to a free escape boundary from a planar source 

in 6 the Mn n 3 on ®' dln,en8ional medium, which is the configuration assumed 
diff^sine 1 reHr° bl f n ’ iS treated under the further assumptions that the 
guiding field a?d It “*" etic turbulence superimposed upon a uniform 
g ng field and that the particles are cosmic rays whose pitch-angles 
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relative to the guiding field undergo scattering that is described by 
the Fokker-Planck formalism. Within this framework, the algebra needed 
to obtain analytic solutions in terms of eigenfunctions becomes 
manageable, because temporal variations are absent. These solutions, 
which embody non-diffusive effects, provide a convenient and well 
understood standard for judging the accuracy with which numerical 
methods describe these effects. 

Some basic aspects of this picture are outlined within these 
proceedings in Paper SH4.1— 4, and this paper deals with the same 
transport equation, invokes the same definitions and notation, and 
insofar as possible, uses the same numbers. 

2. The Method of Eigenfunctions . This approach follows exactly the 
classical treatment in terms of spherical harmonics (Case and Zweifel, 
1967, sec 8.4), except that the distribution function is expanded in 
terms of the scattering eigenfunctions tabulated by Bieber (1977) rather 
than in terms of Legendre polynomials. In the diffusion approximation, 
the distribution function f{z,p} is approximated essentially by the 
first two eigenfunctions. In this case, the isotropic density is finite 
at the surface, and if it is extrapolated beyond the boundary along the 
same uniform slope that applies to diffusion deep within the medium, it 
reaches zero at a distance of (2/3) X outside the boundary. This result 
is exactly the one obtained from the classical Pi approximation. To 
obtain a better description, assume that solutions of the form 

T{y} exp {-z/ A} 

can be expanded as a finite sum of four or more scattering 
eigenfunctions. Then the matrix form of the transport equation can be 
solved only for certain specific values of A that correspond to a 
distinct set of transient eigenfunctions. (Note that the term 
"transient", which derives from classical transport theory, has nothing 
to do with temporal variations.) These are summed in such a way that 
the returning intensity at the surface is as near as possible to zero 
over the outward facing hemisphere. In practice, only the first 
transient is significant. For the examples computed with q = 1.8, which 
are discussed below, the first characteristic length is A = A/ 40. 

3. Numerical Solutions of a Finite-Difference Equation . These 
computations invoke the same formulation and scattering operator as in 
Paper SH4.1-4, except that steady injection was accomplished by 
specifying at a fixed value of z a constant angular distribution. 
Starting from the diffusion solution, computations were carried out for 
a time sufficient to allow a particle moving in a straight line to 
traverse ten times the distance between injection point and boundary. 
At this time, the flux was essentially independent of z, which means 
that the solution was very close to the steady state. (See Earl, 1974b, 
eq. 24.) 

In the figures below, results are plotted for two different 
values of the spatial increment, A z = A/21 (solid curves) and Az * A/210 
(dotted curves), which are, respectively, substantially larger and 
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smaller than the characteristic length A. Note that both increments are 
much less than the value (3/4) A derived by Kota et al. (1982) as a 
condition f °r validity of the diffusive picture. The analytical 
solution described above was computed separately, but it cannot be 
plotted, because its deviations from the dotted curves are smaller than 
the dots. In the figure at the left, the curves at the top give the 
isotropic density and those below give ten times the coefficient f 2 of 
the second scattering eigenfunction, which was obtained by expanding 
computed angular distributions as a series of four eigenfunctions. The 
surface density predicted at z = 0 by the solid curve is ~8% larger than 
expected. This discrepancy can be attributed to inaccuracies 8 in the 
description of non-diffusive effects that can be seen in the curves 
below, which describe an even component that appears in the transient 
an sotropy, but not in the diffusion approximation. Here, the dotted 
curve follows exactly the expected exponential dependence characterized 
by A A/40, while the solid curve, whose spatial resolution is too 
C r\ S V° de ® crlb t thls stron 8 gradient, follows instead an exponential 
which has only 40% of the amplitude expected at the surface, and whose 
decay is characterized by Az = A/21. 

In the figure to the right, the intensity at the surface is 
plotted as a function ot pitch-angle cosine p. The eigenfunction 
soiution is again indistinguishable from the dotted curve, whose rms 
deviation from expected values is 2.2%. Except for the point nearest to 
p - 0 the intensity predicted by the dotted curve within the hemisphere 
P > 0 is essentially zero, while in contrast, the solid curve, which 
? t0 an inappropriate choice of Az, exhibits here a small but 
tinite intensity. Because scattering cannot turn back escaping 

partides beyond the boundary, no returning intensity is expected, and 

nf ' 5 Jff ree ° f abse l nce can be taken as a figure of merit for comparison 
of different calculations. v 


4. The Monte Carlo Method 
statistical model in which 
of particle trajectories 


• The Fokker— Planck formalism derives from a 
scattering results from many tiny deflections 
by turbulent fields. With the aid of 
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computers, these random 
deflections can be 0 90 

simulated directly, and the 
history of individual 0 75 
particles can be traced 0 60 
through space and time. 

Two crucial points of this 0 * 5 
approach are to* choose the 0 30 
amplitudes of the 

deflections to be 0 1S 

consistent with the Fokker- 5 00 
Planck coefficient and o 

to adjust time increments 
so that the systematic 
drift in pitch-angle is 
compatible with the random 
spreading. 

Because of the 
limitations on accuracy 
imposed on this method by 
statistical considerations , 
it is very difficult to 
study the non-dif f usive 
effects considered above. 

However, the figures to the 
right show results obtained in the diffusion regime from a Monte Carlo 
simulation of the Milne problem. Plotted above is the isotropic density 
as a function of distance, which shows the expected uniform slope 
downward to a finite density on the boundary at the right. The arrow 
here designates the point at which the angular distribution plotted 
below was obtained. This angular distribution is in good agreement with 
the one expected from diffusion theory, whose points are designated by + 
symbols • 

5 . Conclusions . The non-dif fusive effects discussed above are 

miniscule, but the precision with which they are described numerically 
leaves little doubt that finite-difference calculations can give valid 
results, provided that their parameters are chosen appropriately. In 
the diffusion regime, Monte Carlo methods offer a useful alternative. 
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